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The efficiency of statistical sampling in broad-histogram Monte Carlo simulations can be considerably im- 
proved by optimizing the simulated extended ensemble for fastest equilibration. Here we describe how a recently 
developed feedback algorithm can be generalized to find optimized sampling distributions for the simulation of 
quantum systems in the context of the stochastic series expansion (SSE) when defining an extended ensemble in 
the expansion order. If the chosen update method is efficient, such as non-local updates for systems undergoing 
a second-order phase transition, the optimized ensemble is characterized by a flat histogram in the expansion 
order if a variable-length formulation of the SSE is used. Whenever the update method suffers from slowdown, 
such as at a first-order phase transition, the feedback algorithm shifts weight towards the expansion orders in 
the transition region, resulting in a non-uniform histogram. 

PACS numbers: 02.70.Rr, 75.10.Hk, 64.60.Cn 



I. INTRODUCTION 

Quantum Monte Carlo techniques are widely employed to 
study equilibrium properties of quantum systems that do not 
suffer from the infamous negative-sign problem which renders 
any efficient simulation impossible yj. Over the last decade 
many technical advances including^the adaptation of non-local 
update techniques 1 2, 3, 4,5,6, 7, 8] from their classical coun- 
terparts J^and the implementation of continuous-time algo- 
rithms IBI H [ipi] - which only later have been adapted to clas- 
sical systems fllj - have made quantum Monte Carlo (QMC) 
simulations almost as powerful as classical simulations, de- 
spite their usually more complex formulation and implemen- 
tation. 

Non-local update techniques efficiently overcome the crit- 
ical slowing-down at many second order-phase transitions 
lfT2i [T3I1 . but do not help to overcome the problem of tun- 
neling out of metastable states at first-order phase transitions 
or the sampling of rough energy landscapes. To alleviate 
these sampling problems for systems with competing, but 
well-separated states a number of techniques have been de- 
veloped that all aim at broadening the range of sampled reac- 
tion coordinates such as parallel tempering / replica exchange 
methods iflill or alternative extended ensemble techniques 
which include multicanonical sampling [15], broad-histogram 
sampling [16] and the Wang-Landau algorithm [17]. Simi- 
lar approaches have been formulated also for QMC simula- 
tions I II8I [Mils El. 

While most of the extended ensemble 
methods sample a flat histogram for a given reaction coor- 
dinate, it has been pointed out that flat-histogram sampling 
generally leads to suboptimal scaling behavior [22J. If ap- 
plied to second-order phase transitions this slowing-down can 
be overcome by introducing non-local update techniques ll23ll . 
In the more general case, it has recently been demonstrated 
that even for local update techniques this slowing-down can 
be overcome by optimizing the sampled statistical ensemble 
||24ll . The idea is to identify bottlenecks of the simulation by 
measuring the local diffusivity of the random walk in a reac- 
tion coordinate and to systematically feedback this informa- 



tion to shift resources (and statistical weight) towards those 
regions where the local diffusivity is suppressed. The result- 
ing histograms are non-uniform sampling distributions which 
are typically tailored to the underlying problem. The feed- 
back algorithm has been applied to a range of classical sys- 
tems with rough energy landscapes including frustrated mag- 
nets 1I25I1 . small proteins ll26ll and dense liquids li27ll . and has 
recently been extended to optimize the simulated tempera- 
ture set in parallel tempering / replica exchange simulation 
schemes ll28ll29ll . 

In this manuscript, we discuss how the ensemble optimiza- 
tion technique can be applied to efficiently sample quantum 
systems in the context of a stochastic series expansion (SSE) 
formulation. Our approach extends previous ideas of adapt- 
ing the Wang-Landau algorithm to quantum systems by ob- 
serving that the stochastic series expansion performs a one- 
dimensional random walk in expansion orders, which simi- 
lar to the classical algorithm can be biased using a general- 
ized density of states (in expansion orders) [20]. It has been 
demonstrated that sampling an extended ensemble can signifi- 
cantly improve the sampling at first-order quantum phase tran- 
sitions lfT9ll20ll . Obtaining a QMC estimate for the general- 
ized density of states allows for the direct calculation of the 
free energy, internal energy, entropy and specific heat for a 
wide range of temperatures. Here we apply the feedback al- 
gorithm to study and further improve the sampling efficiency 
of such extended ensemble quantum simulations. 

For thermal second-order phase transitions where typically 
non-local update algorithms are efficient we find that using 
non-local updates and sampling a flat histogram in expansion 
orders is akeady optimal. However, we find a dependency 
on the underlying representation of the operator string in the 
stochastic series expansion which we discuss in detail. We 
show that the commonly employed fixed length vector repre- 
sentation gives inferior results to a variable length list repre- 
sentation. It is only for the latter the optimal ensemble is just 
a flat histogram in the expansion orders. 

For first-order phase transitions, both thermal and quan- 
tum, we find a significant reweighing when applying the feed- 
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back technique resulting in model-specific histograms. We 
show that these non-uniform sampling distributions further 
improve the simulations when compared to the flat-histogram 
approach, and result in an unbiased distribution of statistical 
errors. 

This paper is organized as follows: In section Ull we out- 
line the stochastic series expansion algorithm and discuss both 
fixed length and variable length representations for the sam- 
pled operator string. We introduce modifications to this algo- 
rithm to sample extended statistical ensembles in section Hill 
and explain the adaptive feedback algorithm that finds the op- 
timized broad-histogram ensemble. Finally, we discuss appli- 
cation of these sampling techniques to a number of quantum 
spin systems in section|lV]and give systematic comparisons to 
flat-histogram sampling methods. We close with a summary 
of our results and in an appendix give a detailed analysis of 
the variable length representation for the operator string and 
the respective update algorithms. 



II. STOCHASTIC SERIES EXPANSIONS 

The fundamental problem for Monte Carlo simulations of 
quantum systems is that the partition function is not a simple 
sum over classical configurations but an operator expression 



Z = Tr cicp{-f3H) , 



(1) 



where H is the Hamiltonian operator, f3 = 1/T the inverse 
temperature (ks — 1), and the trace Tr goes over all states 
in the Hilbert space. The first step of any QMC algorithm 
is thus the mapping of the quantum system to an equivalent 
classical system, and then sample configurations of the con- 
structed classical systems, e.g. a system of world lines. Over 
the years various methods have been developed for this map- 



ping including discrete time 113011 and continuous time ||4 ] pat h 
integrals or the stochastic series expansion (SSE) IBU 13211 . 
While our algorithms can be applied to any of these represen- 
tations 13 3i,l34ll. w hen combined with an extended ensemble 
approach jlSl Il9l l20il . we will concentrate on the stochastic 
series expansion in the following. 

Most commonly, stochastic series expansion is formulated 
in terms of a high-temperature series expansion of the parti- 
tion function 



Z = Tr exp(-/3i/ ) 



/3" 



ri=0 



(2) 

where the expansion coefficients g{n) present a general- 
ized density of states in the expansion order n. This one- 
dimensional representation in expansion orders is similar to 
the classical form of the partition function 



(3) 



and, in fact, we can identify low expansion orders of the ther- 
mal representation (|2|i with high temperature/energy physics 
and higher expansion orders describe low temperature/energy 



physics. In any simulation the expansion has to be truncated 
at some order A which is equivalent to setting a lower temper- 
ature bound. For canonical simulations one typically chooses 
A > 0{N j3) such that contributions from orders n > A are 
negligibly small, and orders n > A are never reached in the 
simulation. 

We can generate a one-dimensional representation of the 
form in Eq. ^ by direct expansion of the Hamiltonian H. 
To this end, we decompose the Hamiltonian into a sum of 
iVb diagonal or off-diagonal non-branching bond terms H — 
X^r*^ Hbi lHH. Inserting this decomposition along with a 
complete set of basis states {|a)} to perform the trace, the 
expression for the partition function becomes 

oo „„ n 

n=0 {q} {S„} ' P=l 

where \a{p)) = 11^=1 |a) denotes the propagated state 
after the action of the first p bond operators on the initial state 
\a) {— |a(0)) = \a{n))). The operator string £"„ is an or- 
dered sequence of n bond-operators Hf,^ , with repetitions al- 
lowed. In order to stochastically sample the terms in the bond- 
operator expansion ^ by their relative weight, we need to in- 
terpret them as probabilities which requires all bond operators 
to be positive semi-definite. For diagonal operators, we can 
simply achieve this by adding an energy offset to each bond 
while for off-diagonal operators there is no general remedy for 
negative weights which is the famous sign-problem for QMC 
calculations. 



A. Fixed length representation 

Most implementations of the SSE algorithm sample the var- 
ious terms in the bond-operator expansion (|4|i by updating the 
operator string 5„ using a fixed length representation. To this 
end, a vector of length A is created which contains n bond 
operators and is filled with A — n identity operators Idp. To 
compensate for the various ways to place the identity opera- 
tors between the bond operators we have to take into account 
a combinatorial factor (^) such that Eq. (|4| becomes 



n=0 {a} {S„} 



^"(A-n)! 
A! 



\{{a{p)\Sr.{p)\a{p~l)) , 
(5) 

where Sn{p) is either a bond operator —Hi,^ of the Hamilto- 
nian decomposition or an identity operator Idp. 

Sampling schemes in SSE algorithms usually consist of two 
distinct types of updates. In a first step, attempts are made to 
change the expansion order n by converting diagonal bond 
operators to identities and vice versa. For the fixed length 
representation, this can be achieved by traversing the operator 
string Sn in ascending order and at each position p propose 
the update if the operator Sn (p) is either a diagonal bond or 
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a) insert operator 

H; H; 
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order n+1 
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b) remove operator 

H; Hi, H: 



L 



order n+1 
order n 



Hi H: 



FIG. 1: Variable length representation of the operator string in a 
stochastic series expansion: Bond operators from the Hamiltonian 
decomposition are marked by filled squares. Between each pair of 
bond operators an empty square denotes the "gap" which in an up- 
date move can be filled with an additional diagonal bond operator as 
shown in the insert move in the top panel. In the reverse move shown 
in the lower panel a diagonal bond operator is removed from the op- 
erator string. The arrow indicates for which position of the operator 
string a subsequent update will be proposed. 



identity operator. The acceptance rates are 



Pace (Idp Hb 
Pace [Hb, 



_ Ndfi {a{p)\Hb^\a{p)) 
A - n 

Id ^ = + 1 

NaPHp)\HbMP)) 



(6) 



where an acceptance probability pacc > 1 is interpreted as 
Pace = 1 as in the usual Metropolis scheme. In both expres- 
sion n denotes the current number of non-identity operators in 
the operator string before the update is accepted/rejected, and 
Nd the number of diagonal bond operators in the Hamiltonian 
decomposition. 

The second part of the update cycle then attempts to switch 
diagonal and off-diagonal bond operators without changing 
the expansion order n. Typically, non-local update techniques 
are employed for this step such as loop updates jHJal> the 
operator loop update tSJ or directed loop updates JTllal- 



B. Variable length representation 

Alternatively to the fixed length representation outlined 
above we can keep the original representation of an operator 
string with a variable number n of (non-identity) bond op- 
erators. Since this list representation of the operator string 
turns out to be the computationally more efficient representa- 
tion in a broad-histogram simulation, we will discuss in the 
following an efficient sequential update algorithm which in- 
serts or removes diagonal bond operators from the operator 
string for this representation, extending on previous work by 
D. C. Handscomb iQ and A. W. Sandvik IM. 



We start with an operator string as illustrated in Fig. [T] 
where non-identity bond operators from the Hamiltonian 
decomposition are depicted by filled squares. These are 
separated by "gaps" into which additional operators can be 
inserted. Starting from the gap preceding the lowest-order 
bond operator we traverse the operator string in ascending 
order, alternatingly updating "gaps" and operators: 



• When updating a gap we try to insert a randomly chosen 
diagonal bond operator, and accept the insertion with 
probability 



Pinscrt 



Nal3{a{p)\Hb^\a(p)) 



1 



(7) 



where n denotes the size of the operator string before 
the update. 

If the insertion succeeds we continue with attempts to 
insert additional operators after the one just inserted, as 
illustrated in the top panel of Fig. [T] and continue insert- 
ing as long as the insertions are accepted. At the first re- 
jected insertion we stop updating this gap and continue 
updating the following operator. 

• When updating a bond operator we attempt to remove 
the current operator if it is diagonal, with probability 



-'remove 



1, 



NaP {a{p)\Hb^\a{p)) 



(8) 



where again n denotes the size of the operator string 
before the update. 

If the removal succeeds we next attempt removing the 
following bond operator and continue removing bond 
operators until an attempt is rejected. We then continue 
updating the following gap. 

The subsequent loop update in the SSE sampling scheme 
which swaps diagonal and off-diagonal bond operators with- 
out changing the expansion order remains unchanged for the 
variable length representation. This algorithm fulfills detailed 
balance as shown in appendix lAl 

Detailed balance is fulfilled not at the level of individual 
update steps, but only at the level of complete sweeps. Hence 
measurements of physical observables have to be done only 
after a full sweep of updates, i.e. a complete traversal of the 
operator string. 



III. OPTIMIZED ENSEMBLES 

The stochastic series expansion is typically performed for a 
canonical ensemble at a fixed inverse temperature j3. From the 
high-temperature expansion of the partition function in Eq. (|2|i 
it is then obvious that the simulated random walk will concen- 
trate on a narrow energy window, or a small number of expan- 
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sion orders characteristic for this energy regime respectively 

p. , oo , 



n=0 



/3 



(9) 

For a system with competing energy scales which are typically 
found in systems undergoing a first-order transition or systems 
with rough energy landscapes the canonical sampling can re- 
sult in dramatic slowing-downs with the random walker stuck 
in one of the energy levels and the tunneling to another en- 
ergy level suppressed by some intermediate energy barrier. To 
overcome this problem parallel tempering or broad-histogram 
algorithms aim at broadening the sampled energy space by in- 
troducing either multiple replicas of the system spread over 
some temperature range or by introducing an extended en- 
semble that creates an additional bias for the random walker 
to tunnel through the intermediate energy barriers. In the fol- 
lowing we concentrate on the idea of generalizing the sampled 
statistical ensemble in the context of a stochastic series expan- 
sion. 

A first step in this direction has been taken by formulating 
the iterative algorithm by Wang and Landau 1 17] to sample a 
broad range of expansion orders in the SSE algorithm i20tl . 
The goal of the algorithm is to approach a flat histogram en- 
semble that samples all expansion orders (up to some cutoff 
A) equally often. Similar to the classical case such uniform 
sampling is achieved by introducing an additional statistical 
weight z«(n) that is inversely proportional to a generalized 
density of states g{n), that is 



Wflat("-) OC 



1 



Wflat [n) OC 



9{n) 
2n+ 1 



fixed length representation, 
variable length representation, (10) 



gin) 

where the additional factor of 2n + 1 is due to the 2n + l local 
steps in one sweep through the operator string in the variable 
length representation (see Sec. IIIBI l in contrast to the fixed 
length representation where the effort per sweep is constant. 
Accordingly the acceptance rates are multiplied by this ad- 
ditional weight factor and setting the inverse temperature to 
(3 = 1. For the fixed-length representation of the operator 
string the acceptance probabilities in Eq. ^ become 



Pace ij^^p 



Pa 



H, 



Idp) = 



w{n + l) Nd{a{p)\Hb^ \a{p)) 

w{n) A — 71 

w(n - 1) A - n + 1 



w(n) Nd{a{p)\HbJa{p)) 



^11) 



and for the variable-length representation of the operator 
string we obtain 



Pace (Dp Hb^) = 
Pace (-ffbp Dp) = 



w{n + l) ^ Nd{a{p)\Hb^ \a(jp)) 

w{n) n + 1 

w{n — 1) n 



w{n) Na{a{p)\HbJa{p)) 



(12) 



In the following we will further expand on this idea and ask 
which weights w{n) are optimal for a given physical sys- 
tem in the sense that the sampled statistical ensemble gives 



fastest equilibration for all thermodynamic observables. For 
classical systems it has been demonstrated that such optimal 
weights exist and generally lead to a non-uniform histogram 
|24]. Typically, the histogram is reweighed and additional re- 
sources (in terms of attempted updates) are shifted towards 
those values of the reaction coordinate (in which the statistical 
ensemble/weights are defined) where bottlenecks of the simu- 
lation occur, typically in the vicinity of a phase transition. The 
optimal weights are systematically approached in Ref. 24 by 
studying a diffusive current of random walkers between the 
extremal values of the reaction coordinate and feeding back 
measurements of the local diffusivity. 

In complete analogy, we can think of the random sampling 
of operator string configurations during a SSE simulation as a 
one-dimensional random walk in expansion orders. To maxi- 
mize equilibration we want this random walker to perform as 
many round-trips between the two extremal expansion orders 
'^min = and Tiinax = A as possible. To identify such round 
trips we add a label to the random walker that keeps track of 
which of the two extremal orders the random walker has vis- 
ited most recently, e.g. "H-" and "-" denote that the walker 
has visited the lowest expansion order n„iin or the highest ex- 
pansion order nmax most recently. The label is switched only 
once the random walker visits the other extremal order Us- 
ing this label we can record two histograms h+ (n) and /i_ (n) 
where after every single update of the operator string we in- 
crement the histogram corresponding to the current label of 
the walker This allows us to measure for each expansion or- 
der what fraction of walkers /(n) = h+{n) / {h^{n) + h^{n)) 
on average are passing by a given expansion order that have 
visited, for instance, the lowest expansion order last. We can 
then define to leading order a current that characterizes the 
diffusion of the random walker from the lowest to the highest 
expansion order 

df 



j = D{n)h(n) 



dn 



(13) 



where D{n) is the local diffusivity, and h{n) = h^{n) + 
h-{n)) the local histogram of visits of the random walker to 
the expansion order n. In order to speed up equilibration we 
want to maximize this diffusive current by varying the sam- 
pled histogram h{n). Following the arguments in Ref. 24 the 
optimal histogram of visits is found to be inversely propor- 
tional to the square root of the local diffusivity D{n) 



h°P\n) cx l/y/D(nj, 



(14) 



where the local diffusivity can be estimated from the current 
histogram of visits by 

n[n) ■ \aj /an\ 

For the example systems discussed in section |IV] we find - 
similar to simulations with classical systems - that this local 
diffusivity depends only weakly on the sampled statistical en- 
semble. The sampled histogram, however, reflects a dramatic 
reweighing of computational resources towards those expan- 
sion orders where the local diffusivity is suppressed, typically 
those orders which can be identified with a thermal or quan- 
tum phase transition as discussed in section HVl 
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A. Feedback algorithm 

In order to obtain a statistical ensemble with weights 
^opt ^^-j which will produce the optimal broad histogram in 
Eq. ( fT4] i we can apply iterations of the feedback algorithm 
outlined in Ref. 24. The algorithm starts from an arbitrary 
broad-histogram ensemble, typically a flat-histogram ensem- 
ble with approximative weights obtained from a few iterations 
of the quantum Wang-Landau algorithm lEoll . The optimized 
ensemble is then systematically approached by feeding back 
the local diffusivity D{n) estimated for the current set of sta- 
tistical weights. To ensure convergence in subsequent feed- 
back iterations, statistical measurements need to be improved 
which we accomplish by increasing the number of simulated 
update sweeps for each iteration. The algorithm can be out- 
lined as follows: 

• Start with some trial weights w{n), 
suchasw(n) w l/g{n). 

• Repeat 

- Reset the histograms 

h{n) = h+{n) = h-{n) = 0. 

- Simulate the system for A'sw sweeps: 

* When traversing the operator string updates 
are accepted according to Eqs. ( II 11121 ). 

* After each step of the traversal update 
the labeled histograms h+{n) and h^{n) 
and the walker's label. 

* Perform a series of loop updates without 
changing the expansion order. 

* Update the global histogram H{n) 
at the end of the sweep. 

- Calculate the fraction f{n) 



{n)+h— (n) 

- Define new statistical weights as 



win) <— 'w{n) 



1 



h+{n) + h-{n) 



dn 



and recalculate the acceptance probabilities. 
- Increase the number of sweeps Ng^ ^ SiVg^. 

• Stop once the histogram h{n) has converged. 

Note that the local histograms h+{n) and h^{n) are up- 
dated after every attempted move when traversing the operator 
string to update the sequence of diagonal bond operators. 

For measurements we need to additionally measure a global 
histogram H{n) after each full sweep. This global histogram 
H{n) is in particular employed to estimate the expansion co- 
efficients g{n) from 



g{n) cx H(n)/w{n), 



(16) 



in terms of the final statistical weights, and normalized via 
g{0) = dn, the dimensionality of the Hilbert space. 



In a fixed length representation, the expectation values of 
the histograms (prior to normalization) are related as h{n) ~ 
A ■ H{n) and there is thus no difference after normaliza- 
tion. In the variable length representation, the average num- 
ber of operator insertion/removal attempts scales with the ini- 
tial expansion order as 2n + 1, so that the number of visits 
h{n) w (2n+ 1) • H{n). 

Based on the values of the expansion coefficients g{n), ob- 
tained using Eq. ( fTSI l, we can calculate thermodynamic prop- 
erties of the system. In particular, the free energy at an inverse 
temperature f3 is obtained as 

F = -ilnZ = -iln^5(n)/3", (17) 

n 

the energy from Eq. Q as 

^=|E'^5W/3""\ (18) 

n 

and the entropy using S — [E — F)/T. In order to ob- 
tain thermal expectation values of an observable A (such as 
e.g. the magnetization), one in addition records separately the 
values of the observable A at each expansion order n, tak- 
ing measurements after each full Monte Carlo step: Denoting 
by Ai the values of the observable A from its measurement 
performed after the i-\h full Monte Carlo step, and by rit the 
value of the expansion order of the corresponding SSE con- 
figuration, the mean values of the observable A at expansion 
order n is given in terms of the global histogram H{ri) as 



A[n) 



1 



H{n) 



(19) 



where the summation is restricted to those values of i, for 
which rii — n. From these mean values A{n) of the ob- 
servable A for each expansion order, one obtains the thermal 
expectation value at an inverse temperature f3 as 



^TrAeM-m = | ^ A{n) g{n)f3'' 



(20) 



Since only a finite number of expansion orders (up to the cut- 
off A) can be treated within the QMC simulations, there exists 
a limited temperature range, over which both the thermody- 
namic quantities and observables can be calculated, while a 
pronounced runaway has been observed beyond this temper- 
ature regime |20]. Furthermore, statistical errors are reliably 
estimated from repeating the full procedure using independent 
random number sequences. 

In our application of the feedback algorithm to a variety 
of quantum spin Hamiltonians we have found that a few itera- 
tions, typically four, are sufficient to obtain convergence of the 
sampled histograms. The initial trial weights are generated by 
running a few iterations of the quantum Wang-Landau algo- 
rithm. For the purpose of the subsequent ensemble optimiza- 
tion it is not necessary for the initial weights to be well con- 
verged to the flat-histogram ensemble, e.g. w{n) — l/g{n), 
but should simply allow to sample a broad histogram over all 
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expansion orders such that round trips between the extremal 
expansion orders occur at a sufficiently high rate such that 
the sampled histogram H{n) and the derivative of the fraction 
f{n) needed for the feedback can be estimated. 



IV. PERFORMANCE STUDIES 

We next present results from a systematic performance 
analysis of the new algorithm for various quantum lattice 
models, including a comparison to flat-histogram techniques. 
In the following, we label results obtained using a flat- 
histogram ensemble by "flat-fixed" and "flat-var" for the fixed 
and variable length operator string SSE representation, re- 
spectively. Results obtained from the optimized ensemble are 
denoted as "opti-fixed" and "opti-var" for the fixed and vari- 
able length representation, respectively. 



A. The spin- 1/2 Heisenberg chain 

We start by considering the case of a spin- 1/2 Heisenberg 
chain. This allows us to compare the various methods first 
in the absence of any phase transition. In later sections, we 
then focus on the additional issues arising from the presence 
of first- and second-order phase transitions. 

We simulate the spin- 1/2 Heisenberg model on a finite 
chain of = 10 sites with periodic boundary conditions, 
described by the Hamiltonian 



(21) 



where J > 0, and Si denotes a spin- 1/2 degree of freedom at 
site i. For this system, we calculated the exact thermal expan- 
sion coefficients gcxact("-) up to order n = A = 500, based on 
a full numerical diagonalization of the Hamiltonian. We used 
these values in order to perform fixed-weights simulations in a 
flat-histogram ensemble. For the fixed length representation, 
we set the weights w{n) — l/5oxact(^^); and for the variable 
length representation, we use w{n) (2n + l)/gexact('T-)- 
In both cases we ran 16 parallel simulations with indepen- 
dent random number streams in order to assess the statistical 
variations due to the QMC sampling. We furthermore per- 
formed up to 4 feedback iteration batches to obtain converged 
optimized ensembles, where the number of MC steps was in- 
creased by a factor of 2 after each feedback step. Again, this 
procedure was repeated with 16 independent random number 
sequences. In all cases, we employed the multi-cluster SSE 
loop update [f^|l2l, flipping each cluster that results from the 
deterministic operator loop construction with probability 1/2. 

In Fig.|2l the averaged fractions f{n) are shown for the var- 
ious algorithms, and Fig. [3] provides the corresponding local 
histograms. Let us first consider the flat-fixed method, which 
is based on a flat histogram ensemble in the fixed length rep- 
resentation. We find that for this method, the fraction f{n) 
in Fig. 12] exhibits steep descents near both ends, n = and 
n = A, of the n-range. As seen from Eq. (flST l. the increased 
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FIG. 2: Average fraction f{n) of random walkers diffusing from 
lowest order (high temperatures) to highest orders (low temperatures) 
in a stochastic series expansion in the inverse temperature of a spin- 
1/2 Heisenberg chain with 10 sites. Data is shown for the optimized 
ensemble with fixed length operator string (opti-fixed) and variable 
length operator string (opti-var), as well as data for the flat histogram 
ensemble for both operator string representations (flat-fixed and flat- 
var). 
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FIG. 3: Sampled local histograms h{n) for the optimized ensemble 
using a variable length operator string (opti-var) and fixed length op- 
erator string (opti-fixed), as well as the corresponding flat-histogram 
methods (flat-fixed and flat-var), in a stochastic series expansion in 
the inverse temperature of a spin- 1/2 Heisenberg chain with 10 sites. 
The inset shows the flat rescaled global histogram (2n + 1) • H{n) 
for the optimized ensemble with variable length operator string. 



derivative df /dn indicates a suppressed diffusivity of the ran- 
dom walker near both ends of the n-range. This slowing-down 
near n = and rt A is however not related to physi- 
cal properties of the system under consideration. Instead, it 
is due to an inefficient local bond-operator insertion/removal 
dynamics in the fixed length operator string representation. 
This is seen from the acceptance probabilities in Eqs. ( fTTT i: 
for small n ^ A the insertion process is suppressed by a 
large denominator, and for n close to A, the removal process 
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Method 


Tup/ 10*' 


'?~down / 10 


Tround/10® 


/ opti — var 
Tround / ^round 


flat-fixed 


1.21(1) 


1.62(1) 


2.83(2) 


6.86(5) 


opti-fixed 


0.892(3) 


0.880(3) 


1.772(5) 


4.30(1) 


flat-var 


0.202(3) 


0.210(3) 


0.412(3) 


1 


opti-var 


0.207(1) 


0.205(2) 


0.412(3) 


1 



TABLE I: Traversal times of random walkers diffusing from low- 
est order to highest orders in a stochastic series expansion in the 
inverse temperature of a spin- 1/2 Heisenberg chain with 10 sites. 
Data is shown for the optimized ensemble with fixed length operator 
string (opti-fixed) and variable length operator string (opti-var), and 
the flat histogram ensemble for both operator string representations 
(flat-fixed and flat-var) in units of lO'' operator insertion/removal at- 
tempts. The statistical error, estimated by performing 16 independent 
runs in each case is given for the least significant digit. 
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is suppressed by a small numerator. The ensemble optimiza- 
tion tries to compensate for this technical inefficiency by shift- 
ing resources towards both ends of the n-range, as seen from 
the histogram shown in Fig. [3] However, this allocation of 
resources near both ends of the n-range is not efficient, as 
the computational effort for performing updates scales pro- 
portional to the expansion order n. It would thus be desir- 
able to overcome the slowing-down of the fixed length oper- 
ator string representation by an improved bond-operator in- 
sertion/removal dynamics. This motivated us to consider the 
variable length operator string SSE representation introduced 
in Sec. II. 

Indeed, changing to the flat-histogram ensembe in the vari- 
able length representation, by using the flat-var method, the 
slowing-down is completely eliminated, as seen from Fig. |2] 
which shows an almost linearly decreasing fraction of the flat- 
var method over the whole n-range. Fig. |2] shows that the 
optimized ensemble method opti-var leads to a very similar 
linear fraction f{n) = 1 — n/A. Indeed, both the flat-var and 
the opti-var method show a flat local histogram h{n), seen in 
in Fig. [3] Interestingly, we thus find that for the Heisenberg 
chain, the flat-var method already provides an optimal ensem- 
ble, in the sense that also the opti-var method leads to a flat 
local histogram h{n). This corresponds in both cases to a uni- 
form diffusivity of the random walker throughout the whole 
n-range. 

Using the variable length representation, the simulation 
of the Heisenberg chain does not suffer anymore from any 
slowing-down, in accordance with the absence of phase 
transitions in this model. The global histogram H{n) w 
/i(7i)/(2n + l) then decreases with n like i7(n) oc l/(2n-|-l). 
This is indeed seen from the inset of Fig. [3] Since the com- 
putational effort of the SSE update methods scales linear with 
n, this implies that an equal amount of resources is devoted to 
each expansion order within the considered n-range. 

The dynamics of the simulation in sweeping through the n- 
range can be quantified in terms of the traversal times of the 
random walker For this purpose, we measure the averaged 
traversal time Tup (rdown) of the random walker to travel from 
the lowest (highest) to the highest (lowest) expansion order, 
in units of single operator insertion/removal attempts. Fur- 
thermore, we denote by Tiound = Tup + Tdown the round-trip 



FIG. 4: Averaged absolute deviation of the calculated expansion co- 
efficients g{n) from the exact result for a spin-1/2 Heisenberg chain 
with 10 sites. Data for the flat histogram ensemble with fixed length 
operator string (flat-fixed) and the optimized ensemble with variable 
length operator string (opti-var) are shown, both averaged over 16 
independent runs. The inset shows the fluctuation around the exact 
result (zero line) of the averaged statistical deviation for the estimate 
calculated in the optimized ensemble. 

time of the random walker In Tab. IIV Al the traversal times 
are compared for the various ensembles and representations, 
averaged over the 16 independent runs in each case. For the 
fixed length representation, we find that in the optimized en- 
semble the random walker's round-trip time Tiound is signif- 
icantly reduced, as compared to the corresponding flat his- 
togram ensemble. Furthermore, for the optimized ensembles. 
Tup and Tdown are similar, meaning that the random walker 
spends about the same number of insertion/removal attempts 
in both directions, whereas in the flat histogram ensembles, 
Tdown is larger than Tup. For the opti-fixed approach we ob- 
tained a larger round-trip time Tiound than for both the flat- 
var and the opti-var method. This shows that the improved 
dynamics of the variable length representation provides a per- 
formance improvement that cannot even be achieved by the 
optimized ensemble using the fixed length approach. That the 
flat-var method is already close to optimal, is reflected by the 
equal traversal times for the flat-var and the opti-var method 
inTab.UVA] 

We next assess the extend to which the improvements in 
the variable length representation reduces the statistical error 
of the expansion coefficients g{n), which are estimated using 
Eq. ( fTSI l from the global histogram H{n) and the employed 
weights w{n), normalized via g{0) = dn, the dimension of 
the Hilbert space (here, dH = 2^=). To quantify the error, we 
calculated the logarithmic deviation 

(51n[g(n)] = In [g(n)] -In [gea;(n)] =\n[g{n)/gex{n)] (22) 

to the exact expansion coefficients, after running each sim- 
ulation for the same, fixed CPU time. In Fig. |4]the absolute 
deviation, averaged over 16 independent runs, (|(51n[g(n)]|) is 
shown vs. n for the opti-var (corresponding in the current case 
to the flat-var method) and the flat-fixed method. We obtain a 
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B. Second-order phase transition 

In our analysis of the spin- 1/2 Heisenberg chain, we found 
that the variable length representation eliminates the techni- 
cal problems with the fixed length representation. Already 
the flat-var approach produces an optimal ensemble with no 
features in the diffusivity, which would otherwise indicate a 
slowing-down of the simulation dynamics in certain regions 
of the n-range. This behavior is expected, as the Heisenberg 
chain does not exhibit any finite temperature phase transition, 
and the multi-cluster updates provide an efficient sampling 
scheme. 

We now consider a system that does exhibit a second-order 
thermal phase transition. In particular, we simulate the spin- 
1/2 Heisenberg model 



FIG. 5: Averaged absolute values of the relative deviation of the 
calculate free energy F for the optimized ensemble using a variable 
length operator string (opti-var) from the exact result for a spin- 1/2 
Heisenberg chain with 10 sites. The data is averaged over 16 in- 
dependent runs. The inset shows the averaged relative deviation of 
the estimated free energy from flat histogram simulations using fixed 
length operator strings (flat-fixed), compared to the optimized en- 
semble method (opti-var). 



significant overall improvement using the opti-var ensemble, 
in particular at low values of n. The inset of Fig.lH showing 
the averaged deviation {5h\[g{n)]) for the opti-var method, 
verifies the agreement within statistical errors of the QMC re- 
sults with the exact expansion coefficients. 

We find the averaged absolute deviation to scale 
(|(51n[g(ri)]|) oc -n}/"^ for the opti-var method. This relates to 
the approximate H{n) (x l/n scaling of the global histogram, 
as expected from the standard scaling of the MC error with the 
number of samples: sampling the systems at an expansion or- 
der n a number of H{n) cx l/n times, the statistical error 
in g{n) cx H{n)/w{n) is anticipated to scale proportional to 

From the expansion coefficients one can calculate thermo- 
dynamic quantities of the model down to temperatures T„iin, 
that scale as T„iin oc 1/A with the cutoff A |20]. Here, we 
consider temperatures T down to 0.03 J, below which strong 
deviations result due to the finite n-range ll20ll . Based on the 
fixed CPU-time QMC results for g{n), we calculated the free 
energy F from Eq. ( fTTI ). and its deviation SF — F—F^x to the 
exact free energy F^x, obtained from the gex{n) in Eq. (fTTT i 
instead. Figure |5] shows the averaged absolute values of the 
relative deviation {\5F/Fex\) for the opti-var method. Com- 
pared to the flat-fixed method, the opti-var approach shows 
significantly reduced errors in the free energy over the whole 
temperature range, in particular reducing the deviation by an 
order of magnitude at high temperatures, T > J (see the inset 
ofFig.Ell. 



H^jY, S^•S,, (23) 

<hj> 

on a three-dimensional cubic lattice, with an antiferromag- 
netic nearest neighbor coupling J > and periodic bound- 
ary conditions. Here, we focus on a cube of Ng — 8^ sites. 
In the thermodynamic limit, the model exhibits a second- 
order phase transition at a temperature Tc ~ 0.946 ifisll . 
which corresponds on this finite system to an expansion or- 
der of around Uc ~ 900. Since no exact results for gex{n) 
are accessible for this system, we first performed 8 Wang- 
Landau sampling steps starting from a uniform assignment 
g{n) — l,n = 0, A — 4000, in order to obtain a broad 
histogram ensemble. Then, we performed 4 feedback itera- 
tion batches, where the number of MC steps was increased by 
a factor of 2 after each feedback step. Similar to the previ- 
ous case, we performed 16 independent simulations for sta- 
tistical averaging the final results. Using the final estimates 
of the expansion coefficients g{n), we then performed simu- 
lations in the flat-histogram ensemble for comparison to the 
optimized ensembles. We generally employed multi-cluster 
updates, except for the opti-var method, for which also single 
cluster (s.c.) updates were considered. In each single cluster 
update step, only one of the clusters resulting from the deter- 
ministic loop construction is flipped. This was done in order 
to artificially reduce the efficiency of the loop update scheme, 
and allows us to study the optimized ensemble when using a 
less efficient SSE loop update scheme with increased autocor- 
relations near the critical temperature. 

The resulting traversal times for the various methods are 
given in Tab. IIV Bl We find that within the variable length 
representation the traversal times are significantly reduced, 
as compared to the fixed length representation. Furthermore, 
also in the current case of a second-order phase transition, the 
flat-var approach performs almost optimal, leading to a simi- 
lar round-trip time than the opti-var method, alert with a mild 
difference between r^p and Tdown, which become equal in the 
opti-var approach. 

In the following, we concentrate on the opti-var method, 
and compare results obtained using either multi-cluster or 
single-cluster updates in the SSE loop updates. The corre- 
sponding round-trip times in Tab. IIV Bl already indicate that 
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Method 


Tup/lO^ 


''"down / 10 


T'round/lO'' 


/ opti — var 
Tround / ''"round 


flat-fixed 


10(1) 


7.8(5) 


18(1) 


7.2(5) 


opti-fixed 


5.87(5) 


5.72(5) 


11.59(7) 


4.6(1) 


flat-var 


1.08(2) 


1.39(3) 


2.46(2) 


1 


opti-var 


1.252(5) 


1.246(4) 


2.47(6) 


1 


opti-var s.c. 


3.13(5) 


3.12(5) 


6.25(5) 


2.5(1) 



TABLE II: Traversal times of random walkers diffusing from lowest 
order to highest orders in a stochastic series expansion in the inverse 
temperature of a spin- 1/2 Heisenberg cube with 8'^ sites. Data is 
shown for the optimized ensemble with fixed length operator string 
(opti-fixed) and variable length operator string (opti-var), as well as 
data for the flat histogram ensemble for both operator string rep- 
resentations (flat-fixed and flat-var) in units of lO'^ operator inser- 
tion/removal attempts. We employed multi-cluster updates, except 
for the opti-var method, for which using only single cluster (s.c.) up- 
dates was considered. The statistical error, estimated by performing 
16 independent runs in each case is given for the least significant 
digit. 
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FIG. 6: Sampled local histograms h{n) for the optimized ensemble 
using single and multi-cluster updates in a variable length represen- 
tation stochastic series expansion of the spin- 1/2 Heisenberg model 
on a cubic lattice with A'^s = 8'' sites. The inset shows the respective 
fractions of random walkers diffusing from lowest to highest expan- 
sion order. For comparison, a linear decrease f{n) = 1 — n/A is 
denoted by a dashed line. 



using single cluster updates, the random walker's dynamics 
suffers from a slowing-down in its dynamics. 

In Fig. |6l we present the histograms for both approaches, 
the inset showing the corresponding fractions. We find that us- 
ing multi-cluster updates, the opti-var method does not suffer 
any significant slowing-down due to the presence of the phase 
transition, as seen from the almost linear decreasing fraction 
/(n). There is only a mild curvature visible; similarly, the 
histogram h{n) does not exhibit pronounced features for the 
multi-cluster approach. Using single-cluster updates is less 
efficient in the critical temperature regime, and the opti-var 
method compensates for this inefficiency by shifting resources 
into the corresponding ?i-range. This illustrates, that the op- 
timized ensemble method in the variable length representa- 



tion, in combination with efficient cluster update schemes, al- 
lows for efficient simulations of a second-order phase transi- 
tion: the cluster-updates eliminate the critical slowing-down, 
so that in the optimized ensemble, a uniform diffusivity can 
be achieved. Furthermore, we find that also in this case of 
a second-order phase transition, the flat histogram approach 
in the variable length representation performs almost optimal, 
showing very similar round-trip times than the optimal ensem- 
ble method. 



C. Thermal first-order phase transition 

In this section, we apply our analysis to a system under- 
going a first-order thermal phase transition. For this purpose, 
we study a system of hardcore bosons on a two-dimensional 
square lattice, with next-nearest neighbor repulsion 1.36.1 . de- 
scribed by the Hamiltonian 

H = -t'^{alaj + a^jUi) + V2 ^ niUk-fi'^ni (24) 

where a| (a^) creates (annihilates) a hardcore boson at site 
i, rii is the density at this site, t the hopping amplitude be- 
tween nearest neighbor sites, V2 > the next-nearest neigh- 
bor repulsion, and ^ the chemical potential. In the following, 
we consider the case of t/V2 = 0.45 for a half-filled lattice 
at /i/V2 — 2. For this parameter values, the model exhibits 
a first-order phase transition at Tc w 0.4 V2 from a high- 
temperature normal fluid phase to a low-temperature smectic 
phase with stripe order |36]. In Ref. |20], the improvement of 
the flat-fixed method over conventional SSE simulations was 
quantified by studying the number of MC steps required to flip 
the directional orientation of the stripes in the smectic phase 
upon crossing Tc- Here, we revisit this model in order to ana- 
lyze the behavior of the opti-var method for a system undergo- 
ing a first-order phase transition. We consider finite systems 
of A^s = L X L lattice sites with periodic boundary condi- 
tions, and scale the cutoff A — 20L^ in order to cover in each 
case a temperature range down to T « O.2V2 < T^. Since 
for this system exact results for g^^ {n) are not accessible, we 
performed a similar procedure as presented in Sec. lIVBl in or- 
der to reach an optimized ensemble and calculate g{n). For 
the non-local part of the QMC updates, we employed the di- 
rected loop method 101, where the number of worms Nw{n) 
constructed at a given expansion order n was kept constant 
during the fixed weights simulations. The appropriate value 
of Nw{n) for each expansion order n was determined during 
the preceding Wang-Landau part of the simulations, such that 
on average 2n bond operators were visited during the non- 
local update part of each MC step 101 ■ 

Figure I2] shows the local histograms h{n) in the optimized 
ensemble method opti-var for different system sizes. In con- 
trast to the flat-var method, we find that the local histograms 
now develop peaks at an expansion order ric ~ 0.55A, in- 
dicative of the first-order phase transition in the correspond- 
ing temperature regime. The opti-var method shifts resources 
into this region of the n-range, where the local diffusivity, c.f. 
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FIG. 7: Optimized local histograms for a stochastic series expansion 
in the inverse temperature for a system of hardcore bosons with a re- 
pulsive next-nearest neighbor interaction on a square lattice of size 
L X L. With increasing system size a peak evolves in the rescaled 
histogram as the system undergoes a weak thermal first-order tran- 
sition corresponding to a characteristic order Uc ~ 0.55A (with a 
maximal expansion order of A = 20L^). 



L 


Ground /lO 


f lat — var t-t r\(i 
■^round / 


6 


0.465(2) 


0.479(7) 


8 


1.848(7) 


2.05(4) 


10 


5.41(2) 


5.60(3) 


12 


12.93(5) 


16.5(2) 



TABLE 111: Traversal times of random walkers diffusing from lowest 
order to highest orders in a stochastic series expansion in the inverse 
temperature of hardcore bosons with a repulsive next-nearest neigh- 
bor interaction on a square lattice of size L x L. Data is shown for 
variable length operator string in the optimized (opti-var) and flat his- 
togram (flat-var) ensemble, in units of 10*^ operator insertion/removal 
attempts. The statistical error, estimated by performing 16 indepen- 
dent runs in each case is specified for the least significant digit. 



Eq. ( fT4l i. shows a suppression that becomes more pronounced 
upon increasing the system size. In contrast to the previous 
cases, we expect such a behavior from the first-order nature 
of the phase transition in the current case: The optimized en- 
semble compensates for the inefficiency of the SSE updates 
to tunnel between the two coexisting phases at this first-order 
phase transition. 

Let us thus compare the performance of the opti-var method 
to the flat-var approach. Tab. IIV CI shows the corresponding 
round-trip times for various system sizes. We find that the 
opti-var method indeed reduces the round-trip time of the ran- 
dom walkers as compared to the flat-var method. The im- 
provement becomes more pronounced with increasing system 
size, as expected from the increasingly sharpening of the lo- 
cal histogram in Fig. [T] Increasing the sampling inside the 
transition region, the opti-var method thus enhances the over- 
all diffusion of the random walker by pushing it towards this 
bottleneck. 
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FIG. 8: Magnetization vs. magnetic field for the spin-1/2 XXZ 
model on a square lattice with 10 x 10 sites at T/J — 0.1. The 
inset shows the magnetization m{nsh) measured at expansion order 
nsh during an stochastic series expansion simulation using the opti- 
mized ensemble in the variable length operator string representation. 



D. Spinflop transition 

Thus far, we performed series expansions in the inverse 
temperature. Finally, we discuss results obtained from per- 
forming a perturbation expansion in a model parameter ll20tl . 
In particular, we consider the spin-1/2 XXZ model 



SfS^ 



(25) 



<»,i> 



on a square lattice with Ng = LxL sites with periodic bound- 
ary conditions, in a finite magnetic field h at an easy-axis 
anisotropy A = 1.5. The model shows a spinflop transition at 
he ~ 1.83 J, where the magnetization m jumps from m = 
for h < he to a finite value of to « 0.12, c.f. Fig. [8] In 
the following, we analyze the behavior of the opti-var method 
as applied to this strongly first-order phase transition. We fix 
the temperature T — O.IJ, and perform an expansion in the 
magnetic field. For this purpose, we decompose h = ho + 6h, 
where ho = 1.5, and introduce separate bond operators for 
the 5h contribution to the magnetic field. 
In the SSE we perform an expansion 



Z = 



E 

n5h = l 



g{nsh)iShy 



(26) 



in the parameter Sh, from which the expansion coefficients 
g{nsh), nsh = 1, ■■■,-^Sh are estimated |20]. In addition, we 
also measure after each full MC update step the total magne- 
tization m = X^i^i '^i of the system, and bin these values 
according to the current expansion order nsh- This way, we 
obtain a magnetization histogram m{nsh)- The inset of Fig. [8] 
shows an example of a magnetization histogram m{nsh) for 
= 0, ...,Asfi = 500. The magnetization histogram is used 
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FIG. 9: Magnetization, optimized local histogram and diffusivity 
of random walkers in a stochastic series expansion for the spin- 1/2 
XXZ model on a square lattice with 4x4 sites. The diffusivity 
shows suppressions for expansion orders where the magnetization 
shows a larger increase. The optimized ensemble compensates for 
the respective slowing-down by shifting additional resources towards 
these expansion orders (see middle panel). 




expansion order Wg^^ /L 



FIG. 10: Sampled local histograms for the optimized ensemble in a 
stochastic series expansion for the spin- 1/2 XXZ model on a square 
lattice with L x L sites. The perturbative expansion at temperature 
T = O.IJ is carried out in terms of the magnetic field Sh. A pro- 
nounced peak in the histogram evolves at a characteristic order of 
the spinflop transition which scales linearly with the linear system 
size L. The inset shows the system size dependence of the round-trip 
time in a linear-log plot. 



to calculate the field dependent magnetization from 

i^)^Z ^ m{riSh)g{nsh){Shr'\ (27) 

such as shown in the main part of Fig. [8] 

Similar to the previous cases, we used a two-step proce- 
dure to obtain an optimized ensemble. We first performed 8 
Wang-Landau sampling steps starting from a uniform assign- 
ment g{nsfi) = 1, nSf^ = 0, Agh — 500, in order to obtain 
a broad-histogram ensemble. Then, we performed 4 feed- 
back iteration batches, where the number of MC steps was 
increased by a factor of 2 after each feedback step. The mag- 
netization histogram is recorded only during the final step. 
Again, we performed 16 independent simulations for statis- 
tical averaging over the final results. 

Let us first consider the results for a small system with 
L = 4 in the opti-var method. In the top panel of Fig. |9] the 
magnetization histogram for this system is shown. It displays 
more structure than the magnetization histogram in Fig.|8]for 
the L = 10 system. In particular, apart from a first increase 
of m{nsh) near nsh ~ 20 (indicated by the solid arrowed 
vertical line), a series of additional increases in m{nsh) are 
found at higher expansion orders (dashed arrowed vertical 
lines). While the increase near ngh ~ 20, followed by a 
mild plateau of to w 0.125 corresponds to the spinflop tran- 
sition, the features at larger ngh relate to the discrete magne- 
tization steps on a finite lattice, which for L = 4 amount to 
Am = l/L^ = 0.0625. Indeed, for T = the magnetization 
process is a series of steps, with jumps of the magnetization 
by Ato. At T = O.IJ, these steps are less pronounced, but 
still visible in Fig. |9] In the plot of the diffusivity D{nsh) 
for the optimized ensemble, the reduced efficiency of the SSE 



algorithm close to magnetization jumps |01 is clearly visible 
(see bottom panel of Fig. |9]l. As seen from the optimized his- 
togram, shown in the middle panel of Fig. |9] the optimized 
ensemble shifts resources into these regions. For the magne- 
tization increase near ngh ~ 20, related to the spinflop tran- 
sition, these features in D{nsh) and H{nsh) are also visible, 
but less dominant. However, upon increasing the system size, 
this slowing-down of the simulation near the spinflop transi- 
tion becomes more pronounced, and dominates the dynamics 
of the random walker 

Figure [To] shows the histogram for the optimized ensemble 
for larger systems. A pronounced peak develops around an 
expansion order, characteristic of the spinflop transition for 
each system size. This feature dominates the optimized his- 
togram, and the additional magnetization steps at larger ex- 
pansion orders become less relevant, as the magnetization in- 
creases smoothly for h> he- Both the position and the height 
of the spinflop peak scale linearly with the linear system size 
L. The width of the peak does not change significantly upon 
increasing L, indicating that the spinflop transition becomes 
increasingly sharp in h, as L increases. 

In Tab. II V Dl we compare the traversal times of the opti-var 
algorithm to the flat histogram methods flat-fixed and flat-var 
for the L ~ 6 system. In the fixed length flat-fixed method, 
strong differences in Tdown ~ 7riip are found, while in both 
variable length approaches (flat-var and opti-var) the random 
walker moves balanced (rdown ~ Tup)- Moreover, in both 
variable length methods, the round-trip time is significantly 
reduced as compared to the flat-fixed method, with the optinal 
ensemble resulting in the largest overall speedup. This im- 
proved simulation dynamics leads to an equivalent reduction 
of the statistical errors in the expansion coefficients g{nsh)- 
We performed in each case 16 fixed CPU time simulations in 
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Method 




"^down / 10 


"nround / 10 


flat-fixed 

flat-var 

opti-var 


3.81(3) 
2.31(5) 
1.95(2) 


21.7(3) 
2.30(3) 
1.92(2) 


25.5(3) 
4.62(3) 
3.87(2) 



TABLE IV: Traversal times of random walkers diffusing from lowest 
order to highest orders in a stochastic series expansion in the mag- 
netic field close to the spinflop transition on a 6 x 6 lattice. Data 
is shown for the flat histogram ensemble with fixed length operator 
string (flat-fixed) and variable length operator string (flat-var), and 
the optimized ensemble for the variable length operator string rep- 
resentation (opti-var), in units of lO'' operator insertion/removal at- 
tempts. The statistical error, estimated by performing 16 independent 
runs in each case is specified for the least significant digit. 

order to assess the standard statistical error in the g{nsh)- The 
results for the L — 6 system are shown in Fig.[TT] While the 
flat-fixed method suffers from a strong increase of the statis- 
tical error over the n-range that relates to the spinflop tran- 
sition, the optimized ensemble method performs more uni- 
formly over the whole n-range. We find a uniform reduction 
of the statistical error in higher n-range; the size of this reduc- 
tion fits well to the square-root of the round-trip time fraction 
shown in Tab. IIV Dl This relation between the statistical error 
and the inverse square-root of the round-tip time was observed 
also in classical systems [2^ . 

Finally, we assess, how the optimal ensemble method scales 
upon increasing the system size. In the inset of Fig. [TO] the 
round-trip time Tround is shown vs. the linear system size L 
for the opti-var method. The observed dependence fits well 
to an exponential increase of Tround with L, indicative of the 
strong first-order nature of the spinflop transition. While the 
optimized ensemble does thus still suffer from an exponential 
performance reduction, it provides a performance improve- 
ment over the flat-histogram method, also in this case. 



0.02 rr 




expansion order 

FIG. 11: Averaged statistical errors of the expansion coefficients 
gijish) calculated in a stochastic series expansion for the spin-1/2 
XXZ model on a square lattice with 6x6 sites. The statistical errors 
are averaged over 16 independent simulations for a fixed CPU time. 
Results are shown for the flat histogram ensemble with fixed length 
operator string (flat-fixed), variable length operator string (flat-var) 
and the optimized ensemble with variable length operator string 
(opti-var). The characteristic expansion order, corresponding to the 
spinflop transition is indicated by the arrow. 



if a variable length formulation of the SSE is used. This flat 
histogram h{n) counting local updates corresponds to a his- 
togram H{n) cx l/(2n + 1) at the level of full sweeps. In 
such computationally "easy" cases it is thus feasible to di- 
rectly incorporate the optimal ensemble weights of Eq. (fTOl ) 
in the quantum Wang-Landau approach, using the variable 
length representation Eq. ( fT2] ). 



V. CONCLUSIONS 

We presented an application of the optimized broad- 
histogram ensemble method to the simulation of quantum sys- 
tems, in order to increase the performance of previously intro- 
duced extended ensemble methods. We found that within the 
stochastic series expansion approach, the fixed length opera- 
tor string representation suffers from an insufficient local sim- 
ulation dynamics at the end of the range for small and large 
expansion orders n, which can be overcome by formulating 
the algorithm using a variable length operator string. We de- 
rived the appropriate update probabilities for this scheme. We 
then adapted a recently developed feedback algorithm to the 
quantum case, and provided a performance analysis of the re- 
sulting algorithm. 

The analysis of the optimized ensemble approach in the 
variable length representation for the Heisenberg chain and 
the second-order thermal phase transition in the Heisenberg 
model on a cubic lattice, suggests that for quantum models 
showing no, or second-order thermal phase transitions, the 
optimal ensemble is characterized by ?l flat histogram h{n) 



In our analysis of the first-order thermal transition of inter- 
acting hard-core bosons and the spinflop transition, we found 
that the optimized ensemble approach shifts the distribution 
of resources towards the transition region, thus allowing for 
a faster equilibration as compared to flat-histogram methods, 
although the scaling with system size remains exponential 
for strongly first order transitions. This suggests that broad- 
histogram ensemble optimization in a single reaction coordi- 
nate might not provide a general remedy of such exponential 
slowing-down. 
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C(i-1) 



FIG. 12: Illustration of the transition C C' by inserting an op- 
erator at position i. In our choice of indexing the operators, after a 
successful insertion all remaining operators shift indices by 2, i.e.. 
Hb^ in C is identical to Ht^^^ ^' fo'^ all fc > i + 1. The ar- 
rows illustrated the positions of the pointers k for intermediate con- 
figurations. All rejected update probabilities cancel except for the 
insertion/removal probabilities of the new operator (marked grey). 



APPENDIX A: PROOF OF DETAILED BALANCE FOR THE 
VARIABLE LENGTH OPERATOR STRING 
REPRESENTATION 



In the following, we prove detailed balance for the vari- 
able length operator string representation update scheme of 
Sec. IIIBI t herebv extending previous work by D. C. Hand- 
scombiSl and A. W. Sandvik iQ. 



In order to obtain the correct equilibrium distribution of 
configurations, it is sufficient to prove that the transition prob- 
abilities from one state to the next satisfy detailed balance: 



w{C) ■ p{C ^ C") = wiC) ■ piC C), 



(Al) 



with 'w{C) the weight of a configuration C and p{C — > C") 
the transition probability from C to C". Here we prove de- 
tailed balance for the variable length operator string represen- 
tation per sweep. The main idea of the proof is that for each 
sequence of operator insertions and removals in a sweep we 
can construct an inverse sequence, and these satisfy detailed 
balance. A configuration C with n operators is defined as a 
sequence of length 2n + 1 consisting of n operators and n+1 
gaps, including one at the beginning and one at the end of the 
string. Configurations during a sweep contain an additional 
"pointer"fc, 1 < fc < 2n + 1, that labels the position at which 
we currently perform insertion and removal operations. While 
the number of operators of the intermediate configurations be- 
tween a sweep may grow and shrink, the pointer k increases 
by zero (removal) or two (insertion) with each step. We label 
such intermediate configurations as C„(fc). 

For the sake of simplicity we first look at the insertion of 
only one operator at position fc during a sweep, as illustrated 
in Fig. [12] Our configuration C containing n operators and 
n + 1 gaps is thus changed into a configuration C" containing 
n + 1 operators and n + 2 gaps, consecutively labeled by an 
index i,l < i < 2n + 1. We label the probability of inserting 
any operator into a gap of an operator string of length n as 
Note that this probability is independent of the position of 
pointer k and depends only on n. With pj^^ (iJt ) the probability 
of inserting operator Hi, and pll^{k) the one for removing the 
operator at position k (even) we obtain 



■ 



■ 



V 



i+2| 
A A 



mm ■ 



C'{i+1) C{i+2) 



C" 



■ 



J g'j+2| 

C"(j+2) 



FIG. 13: The transition from C to C" via an intermediate visit to 
C'. The two newly inserted operators are marked grey and the ar- 
rows denote some of the intermediate configurations encountered in 
the process in forward (filled arrows) and backward (empty arrows) 
direction. 



1 - p"^ 

^ins 



(j-l)/2 

n (1-K™(2j)) (A2) 



21,^ + 1) 



n+l^ 
ins ) 



Nd denotes the total number of bond operators. For comput- 
ing the probability of an operator removal we take configura- 
tion C and traverse it in reverse order. This corresponds to 
reordering the operator string back to front, but since our al- 
gorithm constructs only diagonal (thus commuting) operators, 
this reordering is legal. The transition probability 



2(r, + l) 



PiC ^c)= Y[{i pzt\m ■ (1 - Pit') 



2(n + l) + l- 



(A3) 



{t-l)/2 

n (1 



K™(2j)) 



shows that for only one operator insertion per sweep all the 
rejected insertion and removal attempts cancel and detailed 
balance is fulfilled for our choice of operator insertion and 
removal probabilities, see equations dTjl, (O: 



pjC ^ C) 
p{C' C) 



n+1 



1, 



n+1 



Nal3{a{p)\Ht-^\a{p-l)) 



AC) 
w{C) ■ 

(A4) 

The removal of one operator works analogously. 

In a next step we consider a sweep with two insertion up- 
dates, see Fig. [T3] Again, we start with a configuration C 
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with n operators. We add an operator at position i as before, 
resulting in a configuration C',^j^i{i + 2) of n + 1 operators 
and a pointer at gap i + 2. We then proceed to a position 
j, z + 2 < i < 2n + 1, where we try another insertion. Once 
accepted we arrive at a configuration C''^2(j + 2)- 

Since the transition probabilities between intermediate con- 
figurations depend on the direction of traversal, we denote 
them by arrows ~]> for the forward or for the backward 
direction. We need to prove p(C C")/p{C" — > C) 
'w{C")/w{C). As we have to visit the intermediate states 
+ 2) and C"^2(j + 2) in order to arrive at the final 
state C", we can split up the transition probabilities: 



p{C- 
P{C" 



C") 



C) 



:l?(C'->C;+i(z + 2)) (A5) 

.l?(a;+i(z + 2)^c"), 

+ (A6) 



As before, probabilities cancel when traversing operators and 
gaps in reverse order: 

P (C;'+2(J + 2) ^ C") = P (C" ^ C;V2(J + !))• (A7) 

In order to establish detailed balance, let us perform one in- 
sertion, then visit state C with one more operator as before 
at the right boundary and bounce back to where we were (this 
move cancels), and proceed with the insertion. This will be 
balanced by a removal move, a move to the right and back, 
and another removal move: 



p{c ^ c") _p{c ^ c) "pic ^ + 1)) 



piC'^C) p(C"^C)^(C;+i(z + 2) 
y(C;+i(z + 2)^C") 



C") 



(A8) 



piC" 

'w{C) 

'w{C) 
w{C') 



^(c;+i(z + 2) 



C") 



P{C" 
p(C'- 



p{C" - 
w{C") 



-a 

C") 
■ C) ~p{C' 
w{C") 



^(c;+i(z + i) 



Again, the removal operation is exactly the same. 

We now construct a setup that allows us to treat more than 
one update per sweep. For this we need to take into ac- 
count that there may be more than one sequence of interme- 
diate insertion and removal operations Sj = {Cji, • • • , Cjp} 
that change an initial configuration C via intermediate con- 
figurations in Sj into a final configuration C . We proceed 
by constructing a reverse sequence S' for each S and show 
that most probabilities cancel. The detailed balance condition 



p{C ^ C')/p{C' ^ C) = for each one of these se- 

quences is established in the same manner as above, with a 
"visit to the right boundary" after each successful insertion or 
removal move. 

Again, we start with a configuration C with n operators. 
We add or remove the first operator to arrive at the first inter- 
mediate state of Sj, Cij at position iij as before. We then 
transverse to the right of that operator string by increasing i, 
inserting and removing operators according to Sj and finally 
arriving at a configuration C" in one sweep. The reverse se- 
quence Sj is defined by the intermediate states of Sj in reverse 
order. Analogously to iA5i we obtain 



p{C^C')/p{C' ^C) 



(A9) 



and therefore 



pjc^C) ^ E.pjc^c) 

p{C' ^ C) 



(AlO) 



E.piC'^c) 

J:^w{C')/w{C)p{C' "iC) _ njjC') 



This framework of traversing the operator string to the end and 
^/^back for each operator is only required for the proof. During 
the calculation we just keep a variable length operator string 



w{C) w{C') w{C) 



C (i + 2)) 

"and a pointer k that marks the position at which the next inser- 
tion or removal operation is performed according to the prob- 
abiUties in equations (|7|) , jSj . 
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